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Abstract 

This  study  consisted  of  two  phases.  During  the  first 
phase.  Lumbar  intervertebral  centra  from  healthy,  adolescent 
Rhesus  Monkeys  (Macaca  Mulatta)  were  subjected  to  axial  com¬ 
pression  loads  of  15  and  30  pounds  for  a  period  of  eight 
hours  and  displacement-time  data  was  gathered.  An  axisymmet- 
ric  finite  element  model  was  used  to  analytically  determine 
material  parameters  describing  the  observed  creep  for  each 

v 

applied  load.  A  Three-parameter  Kelvin  Solid  was  used  to 
represent  the  elastic  and  viscoelastic  response.  Parameters 
were  found  by  adjusting  them  until  the  analytical  displace¬ 
ment-time  response  matched  the  experimental  response.  The 
parameters  determined  characterized  the  initial  elastic 
stiffness,  initial  creep  rate,  and  creep  stiffness. 

Analytical  stress  profiles  in  two  horizontal  planes 
through  the  vertebral  centrum  indicated  a  predominance  of 
stress  in  Cortical  bone  and  the  transition  of  stress  from 
the  Trabecular  bone  to  the  Cortical  bone  as  creep  proceeded. 
According  to  the  model,  the  Centrum  behaved  as  though  the 
Cortex  was  acting  like  a  thin  shell  constraining  the  outward 
flow  of  a  viscoelastic  Trabecular  bone  region. 

During  the  second  phase, ' viscoelastic  constants  deter¬ 
mined  for  the  Trabecular  bone  region  were  incorporated  into 
an  overall  model  of  the  intervertebral  joint  minus  the  Arti¬ 
cular  Facet  Joint  and  associated  spinous  processes. 


Viscoelastic  constants  were  then  found  for  the  disk.  How¬ 
ever,  they  proved  to  be  unreasonably  high,  describing  a  very 
stiff,  creep-resistant  disk.  It  was  determined  that  the 
rheological  model  used,  two  Three-Parameter  Kelvin  Solids 
in  series,  was  inadequate  to  allow  one  to  determine  unique 
values  of  these  parameters  in  the  disk. 

Final  results  indicate  that  an  accurate  model  of  the 
intervertebral  joint  must  incorporate  a  viscoelastic  Centrum. 
Creep  behavior  can  no  longer  be  attributed  solely  to  the 


disk. 


A  FINITE  ELEMENT  ANALYSIS  OF  THE 
CREEP  RESPONSE  OF  LUMBAR  INTER¬ 
VERTEBRAL  JOINTS  IN  THE  RHESUS  MONKEY 


I.  Introduction 


1. 1  Purpose 

The  behavior  and  mechanical  properties  of  intervertebral 
joints  are  of  interest  to  researchers  in  several  areas. 

Among  these  are  the  effects  of  sustained  loading,  vibration, 
and  ejection  forces,  as  might  be  experienced  by  aircrews  in 
high-performance  aircraft,  as  well  as  problems  associated 
with  disc  degeneration. 

Intervertebral  joints  have  long  been  known  to  creep 
under  load.  Previously,  the  disk  has  been  treated  as  the 
sole  medium  of  the  creep  (viscoelastic)  behavior.  This  study 
was  undertaken  to  discover  material  constants  which  may 
account  for  such  behavior  experienced  by  the  bony  portion  of 
the  joint  external  to  the  disk  region.  These  values  were 
then  incorporated  into  a  model  of  the  joint  in  order  to  de¬ 
termine  viscoelastic  constants  for  the  disk.  In  both  phases, 
stress  redistributions  resulting  from  creep  behavior  were 
analyzed  for  the  various  internal  regions.  The  results  of 
this  study  will  be  used  to  further  refine  a  biodynamic  model 
of  the  spine. 


1.2  Anatom- 


A  thorough  understanding  of  this  thesis  requires  a  brief 
description  of  applicable  anatomy.  The  spine  will  first  be 
discussed  followed  by  a  more  detailed  description  of  the 
intervertebral  joint. 

The  Rhesus  Monkey  spine  is  a  complex  structure  composed 
of  a  number  of  mobile  vertebrae.  It  is  divided  into  4  reg¬ 
ions,  as  indicated  in  Fig.  1.2-A.  The  concern  of  this  study 
was  with  joints  in  the  Lumbar  region,  specifically  at  the 
L1-L2  level  (indicating  the  top  lumbar  joint) . 

Each  joint  consists  of  two  vertebral  bodies  superior 
and  inferior  to  a  disk  (Fig.  1.2-B);  two  associated  ligaments 
joining  the  vertebral  bodies  with  the  disk,  one  running  ver¬ 
tically  along  the  posterior  edge,  and  one  running  vertically 
along  the  anterior  edge  of  the  joint;  and  an  Articular  Facet 
Joint. 

The  vertebral  body  consists  of  a  centrum  and  a  set  of 
processes  (Fig.  1.2-C)  which  provide  attachment  points  for 
muscles  and  ligaments.  The  spinous  processes  are  attached 
to  the  centrum  by  two  pedicles  and  extend  rearward.  Posterior 
spinous  processes  from  adjacent  vertebral  bodies  form  the 
Articular  Facet  Joint  near  the  edges  of  the  pedicles.  Two 
lateral  processes,  one  on  each  side  of  the  centrum,  extend 
to  the  sides. 

The  centrum  consists  of  spongy  Trabeculae  or  Trabecular 
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Fig.  1.2-A  (Ref  17) .  Rhesus  Monkey  Spine 


Fig.  1.2-B  (Ref  29).  Sagittal  View  of 

Intervertebral  Joint 


Fig.  1.2-C  (Ref  15).  Superior  View  of 

Intervertebral  Joint 


Fig.  1.2-D  (Ref  29).  Portion  of  Intervertebral  Disk 


bone  (Fig.  1.2-B)  encased  circumferentially  in  a  thin  layer 
of  denser  Cortical  bone  (Cortex) .  The  upper  and  lower  sur¬ 
faces  are  composed  of  still  denser  bone  (relative  to  the 
Trabeculae  and  Cortex)  called  the  Bony  End-plate,  which 
serves  as  a  transition  region  between  the  centrum  and  the 
disk. 

The  disk  is  made  up  of  solid  and  fluid  material  in 
three  regions.  In  the  center  of  the  disk  is  a  region  known 
as  the  Nucleus  Pulposus  (Fig.  1.2-D)  comprised  of  a  loose 
network  of  fine  fibrous  strands  in  a  muco-protein  gel.  Sur¬ 
rounding  this  region  is  the  Annulus  Fibrosus,  composed  of 
layers  of  collagen  fibers  imbedded  in  muco-protein  gel. 

These  fibers  are  arranged  in  a  helicoid  manner,  and  fiber 
directions  are  fairly  uniform  within  each  laminate.  In  the 
human  case,  fiber  principal  directions  alternate  between 
adjacent  laminae  at  a  fairly  constant  angle  with  respect  to 
the  disk  midplane  (White  et  al,  Ref  29) .  The  orientation  of 
these  fibers  as  well  as  the  number  of  laminates  is  not  known 
for  the  Rhesus  Monkey. 

Not  shown  in  the  figure  is  the  cartilagenous  end-plate, 
which  is  comprised  of  a  hyaline  cartilage  structure.  It 
serves  as  a  boundary  region  between  a  disk  and  its  adjacent 
vertebral  body. 

1.3  Background 

Many  studies  have  been  performed  in  order  to  gain  an 


understanding  of  the  elastic  and  viscoelastic  response  of 
intervertebral  joints.  Brown  (1957),  Nachemson  (I960), 

Hirsch  (1965) ,  and  Rolander  (1966)  all  conducted  experiments 
on  the  basic  load-deflection  behavior  of  the  unit.  Kulak, 
et  al  (1976)  showed  that  the  human  disk  behaved  like  a  medium- 
length  thick-walled  tube  subjected  to  internal  pressure  (Ref 
18) .  Lin,  et  al  (1978)  concluded  that  human  joints  exhibit 
a  nearly  linear  stress-strain  relationship  (Ref  20) . 

Evans  (1957)  experimentally  examined  the  relative  con¬ 
tributions  of  human  Trabecular  and  Cortical  bone  in  supporting 
axial  loads  and  concluded  that  the  latter  contributed  more. 
Bell  et  al  (1967) ,  however,  noted  that  the  Cortex  was  too 
thin  to  transfer  the  load.  Rockoff,  et  al  (1969)  determined 
that  Cortical  bone  contributes  45  to  70*  of  peak  strength  of 
human  lumbar  vertebrae  (Ref  23) . 

Virgin  (1951)  examined  the  viscoelastic  response  of 
human  disks  and  determined  that  stress  redistributions  were 
related  to  a  fluid  seepage  phenomena  occurring  among  internal 
regions  of  the  joint  (Ref  28) .  Fung  (1960)  showed  that  stress 
in  some  biological  materials  could  be  separated  into  an  elas¬ 
tic  and  a  viscoelastic  part  (Ref  6) .  Kazarian,  et  al  (1978) 
showed  that  the  creep  response  could  be  modelled  by  a  three- 
parameter  Kelvin  Solid  (Ref  17) .  Burns  (1980)  developed  a 
method  for  determining  values  describing  that  model  (Ref  4) . 

Lin  et  jjil  used  a  three-dimensional  model  of  the  joint 
in  axial  compression  to  show  that  posterior  elements  could 


be  ignored  in  such  an  analysis  (Ref  20) .  Tencer,  et  al.  (1982) 
confirmed  that  same  conclusion  for  all  possible  loading  con¬ 
figurations  except  posterior  shear  and  axial  torque  (Ref  27) . 

Belytschko,  et  al  (1974,  Ref  3)  and  Spilker  (1980,  Ref 
24)  independently  employed  axisymmetric  finite  element  models 
to  study  the  elastic  response  of  human  disks.  Hinrichsen 
(1980)  applied  a  similar  approach  to  determine  the  creep 
response  of  the  human  disk  modelled  as  a  three-parameter 
Kelvin  Solid  (Ref  13) .  Allen  (1981)  used  results  from 
Galante  (1967,  Ref  10)  ,  Lin,  et  aJ  (1978,  Ref  19)  ,  and  Kulak, 
et  .al  (1976,  Ref  18)  to  develop  an  inhomogeneous  model  of 
the  human  disk  incorporating  an  orthotropic  Annulus  Fibrosis 
(Ref  2) .  Furlong  (1981)  examined  the  effects  of  degenerated 
disks  on  the  viscoelastic  response  of  Rhesus  Monkey  joints. 

He  assumed  that  the  creep  response  of  the  joint  could  be 
modelled  by  a  homogeneous,  creeping  disk  (Ref  8) . 

1. 4  Assumptions 

In  this  study,  a  disk  and  the  adjoining  vertebral  cen¬ 
trum  will  be  considered  an  "intervertebral  unit".  This  unit 
is  treated  from  observation  as  axisymmetric  about  a  vertical 
centroidal  axis,  as  was  done  by  Furlong  (Ref  8) .  The  verte¬ 
bral  centrum  is  also  assumed  to  have  symmetry  about  a  hori¬ 
zontal  midplane.  Finally,  the  unit  is  considered  to  have 
symmetry  about  the  horizontal  midplane  of  the  disk.  These 
assumptions  will  permit  a  simplified  analysis  requiring  a 


minimi  i  of  computer  resources. 


Any  observed  creep  is  considered  a  quasi-static  pheno¬ 
mena  wherein  inertia  effects  may  be  ignored.  The  magnitude 
of  deformations  are  assumed  small  enough  that  linear  elastic 
response  may  be  assumed  (Ref  15) .  Resultant  internal  stresses 
are  considered  small  enough  that  yielding  does  not  occur; 
hence,  any  plastic  deformation  will  be  ignored.  The  disk  is 
treated  as  nearly  incompressible,  as  expressed  by  other 
authors  (Refs  2,  3,  13,  and  8) .  Biological  seepage  among 
regions  is  ignored.  To  permit  analysis,  each  region  is  con¬ 
sidered  homogeneous  and  isotropic.  Since  the  relative  stiff¬ 
nesses  are  unknown,  the  Bony  end-plate  and  Cortical  bone  are 
initially  assumed  to  have  identical  material  properties; 
later,  the  former  region  will  be  assumed  stiffer  than  the 
latter.  Both  regions  are  considered  elastic.  Since  the 
anatomical  nature  of  the  annulus  fibrosus  is  also  not  pre¬ 
cisely  known,  the  disk  will  be  treated  in  a  11  smeared-out" 
manner  as  one  homogeneous  unit.  Any  viscoelastic  behavior 
is  attributed  to  both  the  Trabecular  bone  and  the  disk  and 
is  considered  linear  for  the  magnitude  of  deformations  experi¬ 
enced. 

1.5  Approach  and  Presentation 

This  thesis  consisted  of  two  phases.  During  the  first 
phase,  vertebral  centra  were  subjected  to  axial  compressive 
loads  of  15  and  30  pounds.  Data  were  gathered  concerning 


axial  displacement  versus  time.  Loads  were  chosen  to  pro¬ 
vide  a  measurable  creep  response,  and  not  necessarily  to 
represent  realistic  anatomical  loads.  Finite  element  tech¬ 
niques  were  used  to  model  the  centrum  and  determine  material 
constants  which  characterized  the  average  experimental  re¬ 
sponse  for  each  load  value. 

During  the  second  phase,  the  viscoelastic  constants 
derived  in  the  first  phase  were  incorporated ' into  a  finite 
element  model  of  the  intervertebral  joint  minus  the  articular 
facet  joint  and  associated  ligaments.  Viscoelastic  constants 
were  then  determined  for  the  disk. 

Starting  with  a  brief  explanation  of  classical  visco¬ 
elastic  theory,  the  following  chapters  describe  some  funda¬ 
mental  concepts  necessary  for  an  understanding  of  this  thesis. 
A  detailed  discussion  of  each  phase  of  analysis  follows. 


The  or 


2.1  Linear  Viscoelasticit 


The  theory  of  viscoelasticity  accounts  for  accumulating 
strains  in  some  materials  when  subjected  to  a  constant  load 
over  a  period  of  time.  It  also  accounts  for  the  relaxation 
of  stresses  experienced  by  those  same  materials  when  subjected 
to  a  constant  strain  over  a  period  of  time. 

From  a  one-dimensional  viewpoint,  Hooke's  Law  relates 
the  stress  and  elastic  strain  experienced  by  a  material  hav¬ 
ing  a  Young's  Modulus,  E.  Such  a  response  may  be  modelled 
by  a  linear  spring  having  a  "stiffness"  of  E  (Fig.  2.1-A). 

Any  strain  produced  by  the  application  of  stress  would  occur 
instantaneously.  The  value  of  that  strain  would  not  change 
unless  the  value  of  the  stress  changed. 

Such  is  not  the  case  in  materials  which  exhibit  visco¬ 
elastic  or  creep  behavior:  strain  may  change,  even  at  con¬ 
stant  stress  values.  The  behavior  of  viscoelastic  materials 


may  be  modelled  by  various  combinations  of  the  linear  spring 
and  dashpot,  which  represents  the  time-dependent  material 
response.  The  simplest  such  model  is  the  Maxwell  Fluid  (Fig. 
2.1-B)  wherein  the  rheological  elements  are  connected  in 
series.  The  total  strain  e  is  the  sum  of  the  elastic  eg  and 
creep  er  strains: 


The  application  of  stress  would  cause  an  instantaneous  elas¬ 
tic  strain  eE  of  the  spring  followed  in  time  by  an  accumula¬ 
tion  of  creep  strain  ec,  which  is  related  to  stress  by 

o  =  n£c  (2-2) 

where  n  is  the  dashpot  or  fluidity  constant.  In  theory,  this 
creep  strain  would  be  infinite;  hence,  this  model  represents 
a  "fluid"  quite  nicely. 

A  more  realistic  model  for  a  solid  material  is  repre¬ 
sented  by  the  Kelvin  Solid  (Fig.  2.1-C),  in  which  the  two 
elements  are  connected  in  parallel.  In  this  model 

e  =  eE  =  ec  (2-3) 


so  that 


=  os  +  oD 


(2-4) 


or 

a  =  Ee  +  ne  (2-5) 

The  application  of  stress  would  cause  strain  only  after  a 
period  of  time,  as  the  dashpot  gradually  permits  strain  to 
occur;  such  strain  would  proceed  until  the  stiffness  of  the 
parallel  spring  will  no  longer  permit  further  deformation. 
Obviously,  even  this  model  is  somewhat  unrealistic  in  that 
it  does  not  permit  any  elastic  (instantaneous)  strain. 

In  order  to  permit  both  elastic  and  viscoelastic  defor¬ 


mation,  the  Three-Parameter  Kelvin  Solid  (Fig.  2.1-D)  is 
used.  It  will  be  incorporated  into  this  thesis.  In  this 


model 


e  =  eE  +  ec  (2-5a) 

so  that 

a  =  EeE  =  qQec  +  q^  (2-6) 

Instantaneous  elastic  deformation  is  performed  by  the  series 
spring  of  stiffness  E.  Creep  strain  rate  ec  is  limited  by 
the  dashpot  having  a  constant  q^.  Ultimate  strain  is  limited 
by  the  parallel  spring  of  stiffness  qQ.  Figure  2.1-E  shows 
a  typical  creep  curve  for  such  a  model.  Rearranging  Eq  (2-6) 
results  in 

dec  ^  (°  q0eC^  (2-7) 

which  shows  that  incremental  creep  strain  dec  depends  on 
total  accumulated  creep  strain  ec  for  any  incremental  time 
step  dt.  Rearranging  Eq  (2-7)  again  results  in 

deC  _  JL  _  %  (2-8) 

dt  qi  ”  qx  eC 

which  indicates  that  q^^  is  the  fluidity  constant  for  the 
dashpot.  This  parameter  is  also  the  slope  of  the  a  versus 
e  curve,  as  seen  in  Fig.  2.1-F. 

Application  of  such  a  model  to  a  multidimensional  state 
of  stress  is  limited  by  its  inability  to  accurately  describe 
the  complex  constitutive  relation  between  stresses  and 


strains.  However,  finite  element  techniques,  discussed  next, 
will  allow  development  of  a  relationship  similar  to  Eq  (2-7) . 

2 . 2  Finite  Element  Method 

2.2.1  General  Concepts .  The  complexity  of  many  problems 
in  mechanics  will  not  permit  exact  closed-form  solutions.  In 
such  cases,  techniques  have  been  used  to  obtain  exact  force- 
displacement  solutions  at  a  finite  number  of  points  and  approx¬ 
imate  solutions  at  other  points  in  the  body.  The  Finite 
Element  Method  is  one  such  technique. 

This  method  presumes  that  the  body  of  interest  can  be 
subdivided  into  elements  within  which  actual  displacement 
functions  may  be  approximated  by  simple  polynomials.  Figure 
2.2-A  shows  one  such  element.  Other  one-,  two-,  and  three 
dimensional  elements  have  also  been  used  successfully.  Each 
element  has  a  set  of  node  points  and  each  node  is  permitted 
certain  rotational  and/or  translational  degrees  of  freedom. 

In  the  case  of  the  triangular  element  shown,  each  node  has 
two  orthogonal  translational  degrees  of  freedom. 

Element  displacements  may  be  written  as  a  polynomial 
function  of  these  nodal  degrees  of  freedom.  The  choice  of 
displacement  function  within  each  element  must  satisfy  three 
requirements.  First,  displacements  must  be  continuous  be¬ 
tween  elements.  Second,  any  rigid  body  motion  of  the  element 
must  not  produce  strain  within  the  element.  Finally,  as  the 
size  of  the  element  is  reduced  to  zero,  the  strain  within  it 
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must  approach  a  constant  value. 

Element  displacements  are  characterized  by  the  total 
set  of  nodal  degrees  of  freedom  {u}.  In  the  case  of  the 
element  shown 

i  T 


{u}  =  {u1  v1  u2  v 2  u3  v 3 } 


(2-9) 

In  this  element,  the  displacement  polynomials  are  written  as 

( 2-10a) 
(2-10b) 


u  =  a^  +  a2r  +  a^z 


v  =  a.  +  arr  +  a,z 
4  5  6 


where  the  a.^  are  constants  to  be  determined. 

For  the  element  shown,  the  nodal  displacements  thus  become 
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( 2-llb) 


Solving  Eq  2-3a  and  Eq  2-3b  for  the  constants  a^  and  substi¬ 
tuting  them  into  Eq  2-?a  and  Eq  2-2b,  respectively,  yields 
an  expression  for  internal  displacements  in  terms  of  nodal 
displacements : 


U  -  Vj 

+  N2U2 

+  N3u3 

( 2-1 2a) 

l 

* 

V  =  NjVl 

+  N2v2 

+  t*3v3 

(2- 12b) 
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where 


Hi  =  2A  <ai  +  bir  +  °iz) 


(2-13) 


a  .  =  r  .  z  -  r  z  . 
i  j  m  m  j 


(2-14) 


b.  =  z  .  -  z 
i  j  m 


(2-15) 


c.  =  r  -  r . 
i  m  j 


(2-16) 


and  the  area  of  the  triangle  is 


1  ri  zi 
A  =  det  1  r2  z2 

1  r3  z3 


(2-17) 


The  indices  i,  j,  and  m  are  permuted  cyclically  in  the  order 
given  and  are  assigned  integer  values  from  1  to  3.  The  pro¬ 
blem  therefore  is  to  compute  a  vector  of  nodal  displacements 
after  which  interval  displacements  may  be  computed. 

2.2.2  Analysis ♦  Problems  of  a  three-dimensional  nature 
in  elasticity  can  be  simplified  by  expressing  the  stress 
tensor  as  the  sum  of  a  hydrostatic  (volumetric)  and  a  devia- 
toric  stress  tensor,  respectively: 


oo  a 
x  xy  xz 


a  „  o  o 
zx  zy  z 


s  o  o 


0  0  0  =  O  S  O  +  0  0  0 

yx  y  yz  yx  y  yz 


o  o  s 


oo  o 
x  xy  xz 


(2-18) 


o  o  o 
zx  zy  z 


In  this  relationship,  s  is  the  hydrostatic  str . ;s  given  by 
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s  =  (ox  +  ay  +  oz)/3  (2-19) 

and 


a  -  s  +  a 
x  x 


a  =  s  +  a 

y  y 


a  =  s  +  a 
z  z 


(2-20) 


t* 


Hydrostatic  and  deviatoric  strain  tensors  are  related  to  the 
total  strain  tensor  in  the  same  manner. 

The  viscoelastic  analysis  used  in  this  thesis  (Ref  6) 
will  assume  that  the  total  strain  at  any  time  can  be  expressed 
as  the  sum  of  hydrostatic  elastic  strain  and  deviatoric  elas¬ 
tic/viscoelastic  strain.  The  total  solution  will  thus  be 
obtained  by  superimposing  an  elastic  and  a  viscoelastic  solu¬ 
tion. 

The  elastic  solution  is  obtained  by  solving 

[K] {u}  =  {P}  +  {Q}  (2-21) 

for  {u},  the  nodal  degrees  of  freedom,  where  [K]  is  the 
elastic  stiffness  matrix,  {P }  is  the  vector  of  applied  nodal 
forces,  and  {Q}  is  the  vector  of  nodal  residual  forces  due 
to  initial  strains.  Appendix  A  describes  the  derivation  of 
[K]  and  {P}  for  an  axisymmetric  problem. 

For  the  elastic/viscoelastic  solution,  the  constitutive 
relationship  for  the  Three-Parameter  Kelvin  Solid  is  inte¬ 
grated  over  time  (Ref  13)  to  give 
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(Aec}  =i^ [! [Al(o)  -  vec>i 


(2-22) 


where  {A£c}  is  a  ve'ctor  of  creep  strains  accumulated  over  a 
small  time  increment  At,  {ec}  is  the  vector  of  total  creep 
strains,  and  [A]  is  given  by 


[A]  = 


1  -vc  _vc  ° 

1  -vc  o 
1  o 


(2-23) 


Sym  2(l+vc) 


Since  the  creep  strain  is  assumed  incompressible,  v  ,  the 

Poisson's  Ratio  for  the  creep  strain  is  0.5. 

Figure  2.2-B  illustrates  the  flow  chart  used  in  the 

elastic/viscoelastic  analysis.  The  quantity  (Aec>  may  be 

considered  an  incremental  initial  strain  {e  }.  Since  the 

o 

initial  strains  enter  into  the  governing  equation  (2-21)  in 
the  residual  force  term  {Q},  Eq  (2-21)  must  be  solved  at  each 
timestep  for  {u}.  Strains  {c}  may  then  be  obtained  from  nodal 
displacements  {u}  by 

U)  =  [B]  {u}  (2-24) 

where 
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Fig.  2.2-B  (Ref  13) .  Flow  Chart  of  Solution  Technique 


As  before,  i,  j,  and  k  are  permuted  cyclically  over  the  in¬ 
tegers  1  to  3,  inclusive. 


Stresses  {0}  are  related  to  the  strains  {e}  by  the  usual 
constitutive  relationship 


{0}  =  [D]  [ { £ }  -  {ec}] 


(2-27) 


where  [D]  is  the  elastic  material  property  matrix  for  an  axi- 
symmetric  body  given  by 


rDi  =  -EU-.v) 

/i  ju\  n 


(l+v) (l-2v) 
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(2-28) 


Finally,  since 

{e}  =  {eE}  +  {ec>  (2-29) 


{ec}  follows  directly.  These  values  of  {0}  and  {ec}  are  then 
substituted  into  Eq  (2-22),  time  is  incremented,  and  the  cycle 
repeats  itself  until  a  maximum  time  is  reached  (Ref  13) . 


Ill .  Centrum  Analysis 

The  first  phase  of  this  study  involved  obtaining  material 
constants  describing  the  creep  response  of  the  Vertebral  Cen¬ 
trum.  These  viscoelastic  characteristics  were  assumed  to  be 
attributable  exclusively  to  the  Trabecular  region. 

During  this  phase,  axial  compression  tests  were  performed 
on  Vertebral  Centra  from  the  L1-L2  level  of  the  spine  of  the 
Rhesus  Monkey  (Macaca  Mulatta) .  Six  specimens  (three  L1-L2 
pairs)  were  loaded  with  15  pounds  (67  Newtons)  and  six  were 
loaded  with  30  pounds  (133  Newtons) .  The  duration  of  loading 
in  both  cases  was  8  hours.  A  16-hour  relaxation  period  fol¬ 
lowed. 

3 . 1  Experimentation 

3.1.1  Equipment.  Figure  3.1-A  shows  the  experimental 
apparatus.  Figure  3.1-B  shows  a  related  block  diagram.  Each 
test  chamber  was  connected  to  a  humidifier.  Data  from  each 
test  chamber  was  passed  first  to  one  of  a  set  of  four  Dana 
DC  amplifiers,  then  to  an  Analog-to-Digital  Convertor  unit, 
and  finally  to  a  Texas  Instruments  Silent  700  ASR  Electronic 
Data  Terminal . 

Each  test  chamber  (Fig.  3.1-C),  incorporated  a  DC-DC 


Linear  Variable  Displacement  Transducer  (LVDT) ,  and  each 
test  cell  was  partially  enclosed  during  testing.  The  front 
of  the  unit  is  shown  open  for  viewing  in  this  figure.  To 
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preserve  the  specimens  during  testing,  100%  humid  air  was 
provided  by  the  commercial  cool-mist  humidifier. 

Voltage  from  the  LVDT  was  amplified  at  a  gain  of  500  by 
an  amplifier.  Sampling  interval  was  adjustable  and  amplified 
signals  were  digitized,  that  is,  converted  to  raw  displace¬ 
ment,  by  the  convertor  unit.  Raw  displacements  were  then 
displayed  and  recorded  on  cassette  tape  by  the  data  terminal. 
Four  data  streams  were  available  and  selectable  at  the  con¬ 
vertor  unit.  Only  two  data  streams  were  used,  however.  In 
this  manner,  two  tests  could  be  run  simultaneously. 

3.1.2  Procedure.  Spines  from  adolescent  Rhesus  Monkeys 
(Macaca  Mulatta)  were  radiographed  to  determine  unacceptable 
abnormalities.  An  L1-L2  unit  was  excised  from  a  selected 
spine  using  a  scalpel  by  cutting  through  the  end  disks  and  by 
disarticulating  the  facet  joints.  Each  unit  was  then  separ¬ 
ated  into  an  Ll  and  an  L2  segment  by  cutting  through  the 
center  disk  midplane  using  the  scalpel.  All  soft  tissue  and 
musculature  were  removed  as  well  as  the  cartilagenous  end- 
plate.  Processes  were  removed  by  cutting  through  the  pedicles 
flush  with  the  centrum  using  a  low-speed  diamond  saw.  The 
prepared  specimen  was  then  frozen  in  a  normal  saline  solution 
to  preserve  it  for  testing. 

Prior  to  testing,  the  specimen  was  thawed  and  then  the 
superior  and  inferior  surfaces  were  photographed  at  a  magni¬ 
fication  power  of  two.  These  photographs  were  used  to  deter¬ 
mine  the  surface  areas,  as  will  be  discussed  in  the  following 
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section.  Then,  the  end  surfaces  of  the  specimen  was  potted 
in  dental  acrylic  using  a  special  mounting  device.  The  pur¬ 
pose  of  the  acrylic  was  to  provide  a  smooth,  flat  surface 
upon  which  the  load  could  be  evenly  distributed.  Figure 
3.1-D  shows  a  specimen  just  after  being  placed  in  one  end  of 
the  device.  Figure  3.1-E  shows  the  final  configuration  with 
the  specimen  in  the  center.  The  potted  specimen  was  removed 
once  the  acrylic  hardened. 

Before  testing  could  begin,  the  LVDT  was  zeroed  and  the 
chamber  displacement  calibrated.  Two  cylinders  were  used  for 
this  task  (Fig.  3.1-F);  one  cylinder  was  1/8  inch  longer  than 
the  other.  With  the  longer  cylinder  in  place  in  the  test 
cell,  the  LVDT  attachment  bolts  were  loosened  and  the  LVDT 
moved  vertically  until  zero  was  displayed  on  a  selected  data 
stream  at  the  terminal.  The  bolts  were  then  retightened. 

Once  zeroing  was  completed,  the  longer  cylinder  was  then 
removed  and  the  shorter  cylinder  pl_aced  in  the  cell.  The  new 
output  was  recorded  after  which  that  cylinder  was  replaced  by 
the  first  cylinder  and  a  third  reading  taken.  In  this  manner, 
outputs  were  compared  to  determine  a  value  representing  a  dis¬ 
placement  of  1/8  inch  as  follows.  The  readings  from  the 
longer  cylinder  (the  one  measured  twice)  were  averaged  and 
the  result  subtracted  from  the  reading  using  the  shorter  cy¬ 
linder  (the  one  measured  only  once) .  This  measured  difference 
thus  represented  1/8  (0.1250)  inch.  All  test  data  was  first 
adjusted  for  zero  and  then  multiplied  by  a  factor  based  on 
this  difference: 
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Calibration  Factor  =  - - -  (3-1) 

|Measured  Dif ference | 10“^ 

The  factor  10“^  was  required  in  Eq  (3-1)  to  convert  this 
measured  difference  to  an  actual  displacement.  The  absolute 
value  is  required  to  convert  negative  differences  to  absolute 
differences.  Appendix  C  presents  an  example  of  this  calibra¬ 
tion  procedure. 

Once  calibration  was  completed,  the  specimen  was  mounted 
in  the  loading  chamber  (Fig.  3.1-G)  and  the  test  commenced. 

The  recorder  on  the  terminal  was  started,  and  the  load  was 
added.  Sampling  times  were  approximately  once  every  second 
for  the  first  ten  minutes,  once  every  ten  seconds  for  the 
following  ten  minutes,  and  then  once  every  minute  for  the 
remainder  of  the  test. 

After  eight  hours,  the  load  was  removed  and  the  above 
sampling  sequence  repeated.  After  another  sixteen  hours, 
the  test  was  concluded  and  the  specimen  removed.  Recorded 
data  was  then  passed  to  the  computer  and  a  displacement- time 
plot  generated  for  the  test  using  an  existing  plotting  pro¬ 
gram. 

Vernier  Calipers  were  used  to  measure  the  specimen  later¬ 
ally  and  longitudinally,  after  which  the  specimen  was  cut 
sagitally  (front  to  rear)  using  a  low-speed  saw  and  a  1  mm 
slice  removed.  This  slice  was  mounted  in  a  dissecting  micro¬ 
scope.  The  Cortex  and  Bony  end-plate  were  then  measured 
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using  a  transparent  graduated  reticle. 


3.1.3  Results .  Appendix  B  contains  the  experimental 
results  for  each  test.  Displacement-time  curves  varied  from 
test  to  test,  especially  between  vertebrae  from  different 
animals.  These  variations  were  due  to  differences  in  size, 
condition,  and  preparation  of  the  specimens. 

Discrete  time  values  were  selected  from  the  first  10,000 
seconds  of  the  loading  period  and  formed  a  basis  for  two 
average  experimental  data  sets,  one  for  each  load  (15  and  30 
pounds).  Only  the  first  10,000  seconds  were  chosen  in  order 
to  preserve  computational  resources  during  the  analytical 
phase  when  these  displacement- time  values  were  matched.  For 
each  test,  displacements  were  extracted  corresponding  to 
these  times.  Then,  for  all  tests  having  a  given  load,  either 
15  or  30  pounds,  these  values  were  averaged  at  each  time  to 
produce  a  set  of  average  experimental  displacement-time  values. 
Two  such  sets  were  thus  produced,  one  for  the  15  pound  load 
and  one  for  the  30  pound  load.  Finally,  each  average  experi¬ 
mental  displacement  was  divided  by  two.  This  was  done  in 
order  that  these  displacement  sets  could  be  compared  to  the 
respective  analytically-generated  sets  of  the  finite  element 
mesh  used,  which  assumed  midplane  symmetry  and  hence  only  dis¬ 
placed  half  of  what  the  full-size  centrum  would  displace. 

Figure  3.1-H  illustrates  these  displacements  for  the  15 
pound  load  case  and  Fig.  3.1-1  shows  the  set  for  the  30  pound 
case.  The  fun j  ;ion 


u ( t)  =  Cx  +  C2t  -  C3e  Xt  (3-2) 

where  u  represents  axial  displacement  in  Cm  and  t  represents 
time  in  seconds,  was  determined  to  provide  the  best  equation 
representing  a  curve  through  these  points.  The  constants 
C^,  C2,  C^,  and  X  for  each  load  case  are  listed  in  Table  3-1. 

TABLE  3-1 

Constants  of  Curve-fit  Equation 
Representing  Model  Displacements 


Constant 

15 

Load 

(lb) 

30 

C1 

1.935  x 

10-2 

2.040  x 

(N 

1 

O 

rH 

C2 

1.900  x 

10- 7 

2.140  x 

10-7 

C3 

9.800  x 

10-3 

5.600  x 

10-3 

-4 

-3 

A 

9.854  x 

10 

1.115  x 

10 

Figures  3. 

1-H  and  3.1-1 

show  that 

most  of  the 

!  creep 

occurred  in  the 

first  2000  seconds  of 

the 

loading 

period  for 

both  loads.  The  initial  slope  of  the  creep  curve  for  the  15 
pound  load  case  was  higher  than  in  the  30  pound  load  case, 
indicating  a  higher  initial  creep  rate.  However,  as  time 
progressed  past  about  4000  seconds,  the  slopes  were  approxi¬ 
mately  equal.  This  fact  would  indicate  that  the  creep  rate 
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depended  on  load  only  initially.  Also,  initial  elastic 
response  for  the  30  pound  case  was  less  than  twice  that  for 
the  15  pound  case,  as  it  would  be  if  the  elastic  stiffness 
was  independent  of  load.  Apparently,  then,  the  initial  elas¬ 
tic  response  was  a  non-linear  phenomena  with  load  as  was  the 
creep  response.  Finally,  the  final  displacement  at  10,000 
seconds  relative  to  the  initial  elastic  displacement  was 
higher  in  the  15  pound  load  case,  also  indicating  that  the 
stiffness  for  the  creep  deformation  depended  on  load. 

3. 2  Analysis 

3.2.1  Development  of  Model.  In  developing  the  finite 
element  mesh,  two  factors  were  considered.  First,  the  exter¬ 
nal  shape  and  dimensions  were  required  to  accurately  reflect 
the  test  specimens.  In  order  to  assure  that  this  require¬ 
ment  was  met,  average  dimensions  were  used.  The  photographs 
of  the  superior  and  inferior  surfaces  were  used  as  a  basis 
for  establishing  an  end-surface  average  radius:  these  areas 
were  averaged  and  the  radius  computed  for  a  circle  of  equi¬ 
valent  area.  The  point  of  maximum  waisting  (where  the  cen¬ 
trum  was  slenderest)  occurred  at  the  midplane.  The  midplane 
radius  was  determined  such  that  the  ratio  of  the  computed 
end-surface  average  radius  to  this  radius  equaled  the  ratio 
of  corresponding  average  values  obtained  from  the  caliper 
measurements.  Finally,  measured  heights  were  averaged  to 
obtain  a  height  for  the  model. 
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Internal  measurements  for  the  Cortical  and  Bony  end- 
plate  were  averaged  in  the  respective  cases.  Each  region 


was  modelled  with  a  uniform  thickness. 

The  second  consideration  was  that  realistic  regional 
material  properties  be  used.  Properties  for  the  Cortical 
and  Bony  end-plate  were  determined  using  a  process  developed 
by  Furlong  (Ref  8) ,  who  formed  the  proportionality: 

EHuman  Centrum  =  EHnman  Cortex  =  EHuman  Trabeculae 
E  g  ( 3— 

Rhesus  Centrum  Rhesus  Cortex  Rhesus  Trabeculae 

This  was  done  because  the  modulus  of  elasticity,  E,  for  the 
Rhesus  Monkey  Cortical  is  not  known.  The  centrum  values  for 
the  human  were  obtained  from  Belytschko,  et  al  (Ref  3)  and 
for  the  Rhesus  Monkey  from  Kazarian  and  Graves  (Ref  15) . 

A  broad  range  of  values  were  listed  for  the  Rhesus  Mon¬ 
key  modulus  of  elasticity  (Ref  15) .  Hence,  it  was  determined 
that  two  values  of  the  modulus  would  be  used  representing  a 
lower  and  upper  bound,  and  that  viscoelastic  constants  would 
be  determined  in  both  cases.  The  values  determined  were 
3,400  Kp/Sq  cm  and  22,000  Kp/Sq  cm,  respectively,  for  the 
Cortex  modulus  of  elasticity.  Since  the  Bony  end-plate  and 
Cortex  were  assumed  homogeneous,  the  Bony  end-plate  was 
assigned  an  identical  modulus  in  each  case.  Poisson's  ratio 
was  0.25  for  both  regions. 
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Figure  3.2-A  illustrates  the  mesh  used  and  its  rela¬ 
tionship  to  the  actual  Intervertebral  Joint.  Internal  regions 
are  also  shown. 

3.2.2  Solution  Technique.  An  existing  finite  element 
program  developed  by  Hinnerichs  (Ref  12)  was  used  to  deter¬ 
mine  the  material  values  for  this  analysis  as  well  as  the 
follow-on  analysis  of  the  entire  joint.  Values  of  E,  qg, 
and  were  adjusted  between  computer  runs  so  that  the 
average  end-surface  displacement-time  response  of  the  finite 
element  model  matched  the  average  experimental  data. 

3. 3  Results 

As  was  discussed  in  the  previous  section,  two  values  of 
the  modulus  of  elasticity  were  used  for  the  Cortex  and  Bony 
end-plate.  These  cases  will  be  referred  to  as  the  "flexible" 
and  "stiff"  cases,  and  correspond  to  the  low-  and  high-modulus 
values  used. 

It  was  experimentally  observed  that  the  Trabecular  region 
reacted  viscoelastically  as  a  function  of  load.  Thus,  for  the 
flexible  case,  finite  element  solutions  were  obtained  to  match 
two  load  levels,  15  and  30  pounds.  Since  the  overall  analysis 
can  be  observed  at  a  given  load,  it  was  assumed  sufficient  to 
study  the  stiff  case  considering  only  the  15  pound  load. 

3.3.1  Flexible  Case.  Figure  3.3-A  shows  the  solution 
for  the  applied  15  pound  load,  and  Fig.  3.3-B  shows  the  solu¬ 
tion  for  the  30  pound  load  case.  Viscoelastic  constants 


A.  Bony  end-plate 

B.  Cortical  bone 

C.  Trabecular  bone 


Fig.  3.2-A.  Centrum  Mesh 
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Fig.  3.3-A.  Finite  Element  Solution,  15  Pound  Load 
Flexible  Case 
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which  produced  a  model-experimental  curve  agreement  for  the 
r  pound  case  were  E  =  160  Kp/Sq  cm,  qQ  =  39  Kp/Sq  cm,  and 
q^  =  77,400  Kp/Sq  cm-sec.  The  values  for  the  30  pound  case 
were  E  =  218  Kp/Sq  cm,  qg  =  176  Kp/Sq  cm,  and  q^  =  281,400 
Kp/Sq  cm-sec. 

Two  planes  (Fig.  3.3-C)  were  chosen  for  analysis  purposes 
in  order  to  study  the  radial  variation  of  stresses  as  well 
as  the  effects  of  creep  on  the  manner  in  which  stresses  were 
redistributed  radially.  These  planes  were  selected  for  a 
better  appreciation  of  the  effects  on  stress  patterns  brought 
about  by  the  proximity  of  the  Bony  end-plate.  The  two  ele¬ 
ments  were  chosen  to  gain  an  even  greater  understanding  of 
these  effects,  as  well  as  others  to  be  discussed  in  the  fol¬ 
lowing  paragraphs.  For  these  elements,  stress  components 
were  studied  as  they  varied  with  time.  It  was  determined 
that  a  combined  analysis  of  stress  in  a  given  plane  as  well 
as  stress  in  a  given  element  in  that  plane  would  provide  the 
greatest  insight  into  the  nature  of  the  stress  distribution 
in  the  Centrum.  Hereafter,  the  top  plane  and  element  in 
Fig.  3.3-C  will  be  referred  to  as  the  "boundary"  set;  the 
lower  plane  and  element  will  be  referred  to  as  the  "midplane" 
set. 

The  original  finite  element  model  assumed  identical 
properties  for  the  Bony  end-plate  and  Cortex.  However, 
the  former  region  is  known  to  be  stiffer,  although  precisely 
how  much  stiffer  is  not  known.  The  results  using  the 


Boundary 

Element 


homogeneous  assumption  for  these  two  regions  was  therefore 
compared  to  the  results  for  an  assumption  of  inhomogeneity, 
wherein  the  Bony  end-plate  was  assigned  a  modulus  of  elas¬ 
ticity  double  that  of  the  Cortex.  Thus,  the  effect  of  in¬ 
creasing  the  Bony  end-plate  stiffness  (relative  to  the 
Cortical)  on  creep  deformation,  stress  distributions,  and 
stress  redistributions  could  be  analyzed  in  order  to  deter¬ 
mine  if  a  more  precise  stiffness  relationship  between  the 
two  regions  is  required  to  adequately  model  the  Centrum. 

Two  cases  were  therefore  analyzed  for  the  applied  15  pound 
load,  a  "homogeneous"  and  an  "inhomogeneous",  depending  on 
whether  the  Bony  end-plate  was  equally  stiff  or  twice  as 
stiff,  respectively,  as  the  Cortex.  The  analysis  incorpor¬ 
ating  the  30  pound  load  will  hereafter  be  referred  to  as 
the  "higher  load"  case  and  will  assume  that  the  Bony  end- 
plate  and  Cortex  are  homogeneous  only. 

In  each  case,  stresses  were  plotted  as  a  function  of 
radius  and  time  for  the  two  planes.  Two  time  values  were 
selected,  0+  and  1100  seconds,  to  represent  stresses  existing 
immediately  after  the  initial  elastic  deformation  and  after 
some  creep  deformation  has  occurred,  respectively.  To  aid 
in  locating  and  comparing  these  plots,  symbols  are  used 
(Fig.  3. 3-D).  These  will  be  annotated  in  the  upper  right- 
hand  corner  of  the  page  having  a  given  plot  for  a  particular 
case.  Plots  are  grouped  by  stress  component  for  comparison 
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Fig.  3. 3-D.  Symbols  Used  for  Quick  Reference 
to  Stress  Plots 
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purposes  among  the  three  cases. 

Figures  3.3-E  through  3.3-P  show  these  plots  in  the 
boundary  plane.  Generally,  stress  patterns  varied  little 
among  the  cases.  The  highest  stress  was  typically  in  the 
Cortex.  (Note:  the  Cortex  position  is  represented  in  each 
plot  by  the  single  point  having  the  largest  radius;  all 
other  radii  are  points  in  the  Trabecular  region) .  The  effect 
of  creep  was  to  transfer  stress  from  the  Trabecular  region 
to  the  Cortex.  Axial  stress  (Figs  3.3-E  through  3.3-G) 
increased  most  notably  in  the  Cortex.  Radial  stress  (Figs 

3.3- H  through  3.3-J)  increased  in  compression  throughout. 

Hoop  stress  (Figs  3.3-K  through  3.3-M)  increased  in  compres¬ 
sion  in  the  Trabeculae  but  transitioned  from  compression  to 
tension  in  the  Cortex.  Shear  stress  (Figs  3.3-N  through 

3.3- P)  decreased  in  the  Trabeculae  while  increasing  in  the 
Cortex.  These  trends  suggest  that,  in  the  boundary  plane, 
the  Centrum  is  behaving  like  a  thin-walled  cylinder  in  which 
the  Cortex  is  a  shell  constraining  the  Trabecular. 

The  effect  of  increasing  the  load  on  the  stress  in  this 
plane  was  most  noticeable  in  the  hoop  stress  component  (Fig. 

3.3- L).  The  maximum  tension  reached  in  the  Cortex  (at  1100 
seconds)  was  actually  less  than  in  the  homogeneous  case, 
even  though  the  stress  at  0+  seconds  was  greater  in  the 
higher  load  case.  Since  the  creep  rate  was  shown  to  be  less 
in  this  case,  the  stress  redistribution  was  also  less. 
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Fig.  3.3-E.  az  Distribution  in  Boundary  Plane,  Homogeneous  Case 
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Fig.  3.3-G.  az  Distribution  in  Boundary  Plane,  Inhomogeneous  Case 


Fig.  3.3-H.  a  Distribution  in  Boundary  Plane,  Homogeneous  Case 


0+  Sec 
1100  Sec 


Plane,  Higher  Load  Case 


Fig.  3.3-J.  or  Distribution  in  Boundary  Plane,  Inhomogeneous  Case 


Fig.  3.3-K.  aQ  Distribution  in  Boundary  Plane,  Homogeneous  Case 
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Fig.  3.3-M.  a0  Distribution  in  Boundary  Plane,  Inhomogeneous  Case 
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Comparing  the  inhomogeneous  with  the  homogeneous  case 
shows  much  less  effect  on  the  stress  in  this  plane.  For  the 
inhomogeneous  case,  all  Cortex  stress  components  were  slightly 
higher  than  the  values  in  the  homogeneous  case,  except  for 
the  shear  stress.  Apparently,  the  stiff er  Bony  end-plate 
had  a  role  in  transferring  relatively  more  stress  to  the 
Cortex  as  it  deformed. 

Figures  3.3-Q  through  3.3-X  show  the  effects  of  homo¬ 
geneity  and  load  on  the  time  variation  of  stresses  in  the 
boundary  element.  The  manner  in  which  each  stress  component 
varied  with  time  was  generally  the  same  in  all  three  cases. 

The  greatest  effect  was  caused  by  the  higher  load  of  30 
pounds,  which  nearly  doubled  each  stress  component  at  each 
time  value  as  compared  to  the  homogeneous  case. 

Again,  the  inhomogeneous  case  differed  much  less  than 
the  higher  load  case  in  terms  of  its  effect  on  these  time 
response  characteristics.  Generally,  the  stiffer  Bony  end- 
plat.e  caused  slightly  higher  stresses  in  the  element,  paral¬ 
leling  the  results  for  the  entire  boundary  plane. 

Stresses  in  the  midplane,  also  generally  followed  a  con¬ 
sistent  pattern.  Axial  stress  (Fig.  3.3-Y  through  3.3-AA) 
generally  followed  the  radial  variation  pattern  determined 
in  the  boundary  plane,  but  the  variation  of  radial  stress 
(Figs  3.3-AB  through  3.3-AD)  was  exactly  opposite  the  trend 
in  the  boundary  plane:  the  lowest  stress  was  in  the  Cortex. 
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Fig.  3.3-Z.  o  Distribution  in  Midplane,  Higher  Load  Case 


Fig.  3.3-AA.  o„  Distribution  in  Midplane,  Inhomogeneous  Case 


Fig.  3.3-AC.  ar  Distribution  in  Midplane,  Higher  Load  Case 


Fig.  3.3-AD.  crr  Distribution  in  Midplane,  Inhomogeneous  Case 


Since  the  Cortex  is  relatively  free  to  expand  outwardly  in 
this  plane,  radial  stress  applied  by  the  Trabeculae  could 
be  relieved.  This  factor  also  explains  why  the  Cortex  exper¬ 
ienced  a  greater  proportion  of  hoop  stress  in  this  plane 
(Figs  3 . 3-AE  through  3.3-AG),  since  radial  displacement  of 
the  Cortex  would  tend  to  increase  this  stress.  Finally  shear 
stress  (Figs  3.3-AH  through  3.3-AJ)  decreased  with  time  in 
the  Cortex  and  increased  in  the  Trabeculae  afe  a  result  of 
creep,  again  opposite  to  the  trend  in  the  boundary  plane. 

These  differences  in  the  stress  variations  in  comparing 
the  two  planes  suggest  that  the  Bony  end-plate  has  a  dual 
role  in  helping  to  restrain  the  radial  displacement  of  the 
Cortex,  as  well  as  transferring  moment  to  it,  as  seen  by  the 
rise  in  local  axial,  radial,  and  shear  stresses.  This  moment 
effect  was  even  more  pronounced  when  the  Bony  end-plate  was 
made  stiffer.  For  planes  nearer  the  midplane  of  the  Centrum, 
this  restraining  mechanism  would  not  aid  the  Cortex  allowing 
it  to  bulge  outward  (Fig.  3.3-AK)  and  hence  experience 
greater  hoop  and  radial  stresses  while  relieving  stress  in 
the  Trabeculae.  If  one  compares  the  variation  of  stress  in 
this  plane,  considering  the  higher  load  case  with  the  homo¬ 
geneous  case,  one  notes  that  the  hoop  stress  in  the  Cortex 
at  1100  seconds  for  the  higher  load  (Fig.  3.3-AF)  was  about 
equal  to  the  value  at  that  location  and  time  for  the  homogeneous 
case.  Apparently,  the  hoop  stresses  within  this  plane,  after 
enough  creep  deformation  has  taken  place,  are  independent  of 
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Fig.  3.3-AF.  oQ  Distribution  in  Midplane,  Higher  Load  Case 


Fig.  3.3-AG.  Og  Distribution  in  Midplane,  Inhomogeneous  Case 
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Doubling  the  Bony  end-plate  stiffness  affected  the 
stresses  more  in  the  midplane  that  in  the  boundary  plane. 

Most  apparent  was  the  relatively  small  amount  of  hoop  stress 
redistributed  in  the  Cortex  (Fig.  3.3-AG)  over  the  1100  sec¬ 
ond  period  for  the  inhomogeneous  case.  This  shows  that  the 
effect  of  a  stiffer  Bony  end-plate  has  an  effect  on  stresses 
in  the  far  field  as  well  as  in  the  near  field. 

An  analysis  of  the  stresses  versus  time  for  the  midplane 
element  confirms  this  conclusion,  as  seen  in  Figs  3.3-AL 
through  3.3-AS,  especially  for  the  axial  and  shear  components. 
Over  time,  for  all  three  cases,  the  element  experienced  a 
reduction  in  axial  and  shear  stresses,  while  the  radial  and 
tangential  (hoop)  components  increased.  Most  interesting  is 
the  fact  that  the  radial  and  tangential  components  were  less 
in  this  element  for  the  higher  load  case  than  for  either  15 
pound  load  cases.  As  before,  this  can  be  traced  to  the  radial 
displacement,  as  the  element  lengthened  radially  nearly  10% 
more  under  the  30  pound  load  (Fig.  3.3-AT).  In  this  manner, 
the  higher  axial  load  caused  the  element  to  be  relieved  of 
radial  and  hoop  stress  by  these  radial  displacements,  as  the 
element  achieved  a  state  of  equilibrium. 

The  overall  results  for  the  flexible  analysis  imply 
that  the  effect  of  load  was  to  cause  a  greater  radial  dis¬ 
placement  of  the  Centrum  midplane.  The  assumption  that  the 
Bony  end-plate  and  the  Cortex  are  homogeneous  is  reasonable, 
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since  doubling  the  stiffness  of  the  former  had  little  effect 
on  the  overall  creep  response  of  the  Centrum.  Unless  the 
Bony  end-plate  is  significantly  stiffer  than  was  assumed, 
the  homogeneous  assumption  is  adequate. 

3.3.2  Stiff  Case.  To  account  for  the  upper  bound  on 
the  range  of  values,  the  Cortex  and  Bony  end-plate  were  next 
assigned  an  E  value  much  greater  than  before,  and  a  new  set 
of  viscoelastic  parameters  found  for  the  Trabecular  region. 
These  new  values  were  E  =  98.5  Kp/Sq  cm,  q0  =  14  Kp/Sq  cm, 
and  q-^  =  28,140  Kp/Sq  cm-sec.  Each  value  was  less  than  its 
respective  counterpart  determined  using  the  flexible  Bony 
end-plate  and  Cortex,  again  underscoring  the  influence  of 
the  two  regions  on  the  creep  response  of  the  Centrum. 

To  provide  a  check  on  the  results,  Eq  (3-3)  of  section 
3.2.1  was  used  to  compare  the  known  ratio  of  Human  to  Rhesus 
Monkey  Centrum  moduli  of  elasticity  to  the  analytically 
determined  ratio  of  Human  to  Rhesus  Monkey  Trabeculae  moduli 
of  elasticity.  Using  the  information  in  that  section  as  well 
as  a  Human  Trabeculae  E  value  obtained  from  Belytschko,  et  al. 
(Ref  3),  the  two  ratios  were  6.2  for  the  Centra  and  7.6  for 
the  Trabecular  regions.  Hence,  Eq  (3-3)  is  a  useful  method 
to  approximate  Rhesus  Monkey  moduli  of  elasticity,  in  the 
absence  of  precise  values. 

Figure  3.3-AU  shows  the  analytical  solution  for  average 
top-surface  displacement.  To  avoid  repetition  of  the  analyses 
of  the  flexible  case,  only  the  15  pound  load  data  was  matched. 
Figures  3.3-AV  through  3.3-BC  show  the  stress  distributions 
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Fig.  3.3-AV.  ar  Distribution  in  Boundary  Plane,  Stiff  Case 


Fig.  3.3-AW.  or  Distribution  in  Midplane,  Stiff  Case 


Fig.  3.3-AX.  a_  Distribution  in  Boundary  Plane,  Stiff  Case 


Fig.  3.3-AZ.  a0  Distribution  in  Boundary  Plane,  Stiff  Case 
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Fig.  3.3-BA.  Oq  Distribution  in  Midplane,  Stiff  Case 


Fig.  3.3-BC.  arz  Distribution  in  Midplane,  Stiff  Case 


in  the  boundary  plane  and  midplane,  where  the  plots  are 
again  organized  by  stress  component  for  comparison  purposes. 
The  same  general  trends  existed  in  each  plane  regarding  radial 
stress  variation  and  stress  redistribution,  as  existed  in 
those  planes  for  the  homogeneous  case  (flexible  Cortex)  dis¬ 
cussed  in  the  previous  section.  Differences  were  in  the 
relative  contribution  of  Cortex  and  Trabeculae  in  carrying 
loads.  In  the  stiff  case,  the  Cortex  carried  a  proportion¬ 
ately  larger  share  of  radial  and  tangential  stresses  in  the 
boundary  plane  (Figs  3.3-AV  and  3.3-AZ).  Trabeculae  stresses 
also  tended  to  be  more  evenly  distributed  in  both  planes. 

The  midplane  element  experienced  the  same  general  stress¬ 
time  response  (Figs  3.3-FD  through  3.3-BG),  but  the  axial 
stress  for  the  boundary  element  (Fig.  3.3-BE)  decreased  with 
time  instead  of  increasing,  as  it  did  for  the  homogeneous 
case  of  the  previous  section.  The  stiffer  Bony  end-plate/ 
Cortex  combination  allowed  the  axial  stress  to  relax  under¬ 
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IV.  Joint  Analysis 


The  second  phase  of  analysis  for  this  investigation  in¬ 
volved  obtaining  viscoelastic  constants  for  the  disk  modelled 
as  a  homogeneous  region.  These  constants  described  the  aver¬ 
age  creep  response  reported  by  Furlong  (Ref  8) ,  who  obtained 
creep  data  on  healthy  intervertebral  units  from  the  L1-L2 
level  of  the  Rhesus  Monkey  spine.  He  conducted  his  tests  by 
subjecting  specimens  to  an  axial  compressive  load  of  15 
pounds  (67  Newtons)  for  8  hours  followed  by  a  16-hour  relaxa¬ 
tion  period.  Specimens  were  obtained  by  cutting  the  joints 
through  the  disk  midplane  and  through  adjacent  Centrum  mid¬ 
planes  (Fig.  4.1-A). 

4. 1  Finite  Element  Model 

The  finite  element  mesh  used  for  this  analysis  was  ob¬ 
tained  from  Furlong  and  is  also  shown  in  Fig.  4.1-A.  The 
material  properties  (E  and  u)  used  for  the  Trabeculae,  Bony 
end-plate,  and  Cortex  were  identical  to  those  used  in  the 
Centrum  analysis  previously  discussed  in  this  thesis.  The 
Bony  end-plate  and  Cortex  were  assumed  homogeneous  with 
E  =  3,400  Kp/Sq  cm  and  u  =  0.25.  The  viscoelastic  constants 
(E,  qQ,  and  q^)  for  the  Trabecular  region  were  those  found 
in  the  Centrum  analysis  for  the  flexible  case  at  an  applied 
load  of  15  pounds.  Poisson's  ratio  for  the  disk  v<as  0.48. 

E,  qQ,  and  q^  for  the  disk  were  initially  set  equal  to  the 
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values  Furlong  obtained  for  the  disk  as  he  modelled  it. 


4.2  Film  Data 

One  major  limitation  of  the  method  of  analysis  used  in 
this  study  was  that  one-dimensional  data  (top  surface  dis¬ 
placement)  was  used  to  approximate  properties  describing  a 
two-dimensional  displacement  phenomena.  A  means  was  required 
for  confirming  the  validity  of  the  analytical  results  ob¬ 
tained  from  this  standpoint. 

A  series  of  time-lapse  films  was  used  to  provide  such 
a  check.  These  films  were  taken  by  Furlong  (Ref  8) .  The 
test  cell  was  filmed  during  each  test;  each  film  included  a 
clock  to  indicate  elapsed  time.  The  diameter  of  the  disk 
midplane  (mid-height  of  the  specimen)  was  measured  using  a 
travelling  microscope.  Three  tests  were  chosen  (one  test 
per  film)  and  five  frames  were  selected  from  each  film  to 
obtain  these  measurements. 

Two  major  difficulties  were  encountered  in  taking  these 
measurements.  First,  the  amorphous  nature  of  the  biological 
material  as  well  as  a  lack  of  contrast  made  measurement 
points  hard  to  see  at  times.  Secondly,  the  accuracy  required 
to  measure  the  very  small  displacements  of  the  disks  between 
successive  frames  very  nearly  approached  the  accuracy  limits 
of  the  microscope.  However,  using  a  technique  where  each 
frame  was  measured  several  times,  an  average  measurement  was 
then  obtained  for  that  frame. 
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Only  the  first  2,000  seconds  of  each  test  film  was  used, 
since  that  was  the  limit  of  the  experimental  displacement 
time  data  modelled  in  this  analysis.  Measurements  were  taken 
at  approximately  400  seconds  into  each  test  and  at  400  second 
intervals  thereafter.  Since  frames  were  exposed  once  every 
25  seconds,  these  intervals  could  only  be  approximated. 

Disk  midplane  radii  were  computed  from  each  diameter 
measured  by  simply  dividing  by  two.  The  difference  in  radii 
between  chronologically  successive  frames  thus  represented 
an  unsealed  displacement  over  the  interval.  These  values 
will  hereafter  be  referred  to  as  measured  interval  displace¬ 
ments,  AUr  .  Each  film  yielded  four  such  values. 

These  measured  interval  displacements  were  adjusted  to 
full-scale  interval  displacements,  AUr,  using  a  scale  factor 
obtained  by  dividing  the  measured  disk  radius  for  the  400- 
second  (earliest)  frame  into  the  disk  radius  used  in  the 
finite  element  model  (0.85  cm).  Thus 

0  •  85 

Scale  Factor  =  Disk  radius  at  400  Sec  (4-1) 

The  actual  interval  displacements  were  therefore  expressed  as 

AU_  =  AUr  •  (Scale  Factor)  (4-2) 

L 

Interval  median  times  were  determined  for  each  time 
interval.  These  times  were  based  on  the  400-second  initial 
measurement  as  well  as  its  corresponding  time  interval.  Four 
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average  median  times,  t  ,  were  then  obtained  from  the  four 

m 

time  intervals  in  each  of  the  three  films. 

Finally,  for  each  time  interval,  the  four  values  of  AUr 
from  the  films  were  averaged  to  obtain  an  average  interval 
displacement,  AUr.  The  results  of  this  procedure  are  listed 
in  Tables  4.2-A  through  4.2-D.  Good  agreement  existed 
between  these  values  and  the  analytical  displacements  re¬ 
ported  by  Furlong  (Ref  8)  . 

4 . 3  Results 

Experimental  displacement- time  values  and  the  analyti¬ 
cally  generated  data  are  shown  in  Fig.  4.3-A.  The  displace¬ 
ment  function  described  in  section  3.1.3  (Eq  3-2)  was  again 
used  to  provide  a  curve-fit  to  the  experimental  data.  Con- 

—  6 

stants  defining  this  curve  were  =  0.0312,  C 2  =  9.333  x  10  , 

C3  =  0.0055,  and  X  =  0.0432. 

Viscoelastic  constants  obtained  for  the  disk  were  E  = 
1,000  Kp/Sq  cm,  qg  =  1,000  Kp/Sq  cm,  and  q^  =  99,000  Kp/Sq 
cm-sec.  It  is  the  author's  opinion,  based  on  previous  experi¬ 
mental  results  reported  by  Furlong  (Ref  8) ,  Allen  (Ref  2) , 
and  Hinnrichsen  (Ref  13) ,  that  these  constants  are  unrealis¬ 
tically  high,  indicating  a  disk  which  is  stiffer  and  creeps 
at  a  slower  rate  than  the  Trabecular  region  in  the  Centrum. 
Efforts  to  match  the  experimental  curve  precisely  resulted 
in  viscoelastic  constants  at  least  three  orders  of  magnitude 
larger  than  those  determined.  One  possible  reason  for  this 
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discrepancy  could  be  that  the  rheological  model,  consisting 
of  two  three-parameter  Kelvin  solids  in  series,  used  to 
represent  the  creep  deformation  of  the  overall  joint  will 
not  permit  determination  of  unique  viscoelastic  constants 
for  the  disk  given  similar  constants  for  the  Centrum.  Per¬ 
haps  a  Three-parameter  Kelvin  Solid  in  series  with  a  Kelvin 
Solid  would  yield  better  results.  Further  investigation  is 
required  to  determine  a  viable  model. 

To  provide  a  check  on  the  creep  deformation  experienced 
by  the  model,  the  film  data  determined  in  the  previous  sec¬ 
tion  for  the  radial  displacement  of  the  disk  midplane  was 
compared  to  the  model  deformation  of  a  corresponding  node 
located  at  the  outer  edge  of  the  disk  midplane.  The  results 
of  this  comparison  are  summarized  in  Table  4.3-A.  One  notes 
immediately  that  the  model  displaced  several  orders  of  mag¬ 
nitude  less  than  the  film  measurements  predicted,  further 
implying  that  the  disk  has  been  modelled  as  being  too  stiff. 

Despite  the  limitations  of  the  analytical  results,  it 
was  determined  that  an  analysis  of  stress  patterns  would 
prove  beneficial  in  describing  the  manner  in  which  stresses 
are  distributed  and  redistributed  in  the  intervertebral  joint. 
Two  planes  were  selected  for  this  analysis  (Fig.  4.3-B). 
Figures  4.3-C  through  4.3-E  present  these  plots  for  the  disk 
plane  and  Figures  4.3-G  through  4.3-J  present  them  for  the 
centrum  plane.  In  the  disk  plane,  radial  stress  (Fig.  4.3-C) 


TABLE  4.3-A 


Comparison  of  Disk  Radial  Displacement 
From  Film  Measurements  and  Analytical  Results 


AUr  (cm) 

Median  Time  (sec)  Films  Analytical 


610 

644 

1019 

1065 

1417 

1465 

1821 


0.00270 


0.00224 


0.00196 


0.00143 


0.0000124 


0.0000225 


0.0000084 


1840 


-0.0000012 


*  Experimental  (Ref  8) 

A  Finite  Element  Solution 
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Fig.  4.3-A.  Finite  Element  Solution, 
Joint  Analysis 
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Fig.  4.3-C.  0  Distribution  in  Disk  Plane 


Fig.  4. 3-D.  a7  Distribution  in  Disk  Plane 


varied  little  with  radius  and  generally  increased  in  magni¬ 
tude  as  creep  occurred.  Axial  and  tangential  stresses  (Figs 

4. 3- D  and  4.3-E)  also  varied  little.  Hoop  stress  was  tensile 
near  the  edge  of  the  disk,  indicating  that  the  disk  is  acting 
like  a  shell.  The  greatest  radial  variation  occurred  in  the 
shear  stress  component  (Fig.  4.3-F)  where,  over  time,  stress 
magnitude  generally  increased  for  most  radii.  The  fact  that 
little  stress  redistributed  in  this  plane  can  again  be  traced 
to  the  high  value  determined,  indicating  a  relatively 
creep-resistant  material.  These  patterns  are  thus  in  dis¬ 
agreement  with  prior  studies  (Refs  2,  8,  and  13)  . 

Stresses  in  the  Centrum  plane,  however,  are  more  real¬ 
istic,  parallelling  the  results  for  the  midplane  in  the  Cen¬ 
trum  analysis  of  this  thesis.  Radial  stress  (Fig.  4.3-G) 
was  compressive  and  increased  in  magnitude  as  creep  occurred. 
Axial  stress  (Fig.  4.3-H)  was  highest  in  the  Cortex  (data 
point  having  the  highest  radius) ;  stress  was  transferred  to 
the  Cortex  as  creep  progressed.  Tangential  stress  (Fig. 

4.3- 1)  was  compressive  in  the  Traceculae  and  tensile  in  the 
Cortex;  stress  magnitudes  increased  at  all  radii.  Finally, 
shear  stress  (Fig.  4.3-J)  varied  considerably,  especially  at 
the  interface  between  the  Trabeculae  and  Cortex,  and  also 
increased  in  magnitude  over  time.  These  patterns  reconfirm 
the  results  of  the  Centrum  analysis  and  again  describe  a 
Centrum  whose  Cortex  is  behaving  like  a  thin  shell  constrain¬ 
ing  the  outward  flow  of  the  Trabeculae. 
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V.  Conclusions  and  Recommendations 


Centra  from  Rhesus  Monkey  intervertebral  joints  were 
subjected  to  axial  compressive  loads  for  eight  hours  and 
observed  to  exhibit  viscoelastic  (creep)  behavior.  Using 
an  axisymmetric  finite  element  model  of  the  Centrum  which 
incorporated  a  viscoelastic  Trabecular  bone,  material  para¬ 
meters  quantifying  this  behavior  were  determined  for  applied 
loads  of  15  and  30  pounds  by  matching  one-dimensional  experi¬ 
mental  displacement  data. 

Significant  stress  redistributions  occurred  in  the  Cen¬ 
trum  as  a  result  of  creep  behavior.  It  was  found  that  the 
Trabecular  bone,  as  it  creeps,  transfers  load  to  the  Cortex. 
Thus,  the  Centrum  can  be  characterized  as  having  a  Cortex 
which  is  acting  like  a  thin  shell  constraining  the  outward 
flow  of  the  Trabecular  bone. 

It  was  determined  that  it  is  reasonable  to  consider  the 
Cortex  and  Bony  end-plate  as  a  homogeneous  unit,  in  the 
absence  of  precise  material  properties  for  the  latter  region. 

Load  level  had  an  effect  on  the  material  parameters 
found  in  the  Centrum.  In  general,  all  three  constants  were 
so  dependent,  as  the  overall  creep  response  was  non-linear 
with  load. 

For  the  follow-on  analysis,  viscoelastic  constants 
determined  for  the  Trabecular  bone  were  applied  to  an  overall 
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model  of  the  intervertebral  joint  (minus  the  articular  facet 
joint  and  associated  ligaments)  with  the  object  of  determin¬ 
ing  similar  constants  for  the  disk  region.  Constants  deter¬ 
mined  were  unreasonably  high,  indicating  an  excessively 
stiff,  creep-resistant  disk.  It  was  determined  that  the  use 
of  two  Three-parameter  Kelvin  Solids  in  series  may  be  an 
inadequate  rheological  model,  since  unique  parameters  for 
the  disk  were  unobtainable  by  the  analytical  method  utilized. 
It  is  recommended  that  other  rheological  models  be  considered 
to  ascertain  a  more  viable  model  of  the  overall  joint  creep 
response.  Viscoelastic  constants  for  the  disk  should  then 
be  determined. 

A  crude  optical  technique  was  used  wherein  time-lapse 
film  images  were  measured  to  determine  disk  radial  displace¬ 
ments,  and  provide  a  check  on  the  two-dimensional  creep 
response  of  the  finite  element  model.  It  is  recommended 
th.it  such  a  technique  be  explored  to  ascertain  ways  of 
improving  related  photographic  and  measurement  techniques. 

The  incorporation  of  a  viscoelastic  Trabecular  bone 
region  is  necessary  to  precisely  model  the  time  dependent 
deformation  of  the  intervertebral  joint.  As  such,  the  disk 
cannot  be  considered  the  sole  medium  of  creep  behavior. 
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APPENDIX  A 


Stiffness  Matrix  and  Nodal  Force 
Vector  For  Axisvmmetric  Problems 

The  derivation  of  the  Stiffness  Matrix  and  Nodal  Force 
Vector  proceeds  as  follows  for  an  axisyirunetric  problem. 

A.l  Stiffness  Matrix  (Ref  31) 

The  general  equation  for  the  stiffness  matrix  of  an  ele 
ment  is  given  by 

[k]  =  /  [b]T[D] [B]d (volume)  (A-l 

V 

For  an  axisyirunetric  problem,  this  integration  must  proceed 
over  the  entire  ring.  Hence, 

[k]  =  2rr/  [B]T[Dj  [B]rdrdz  (A-2! 

Summing  the  contributions  from  each  of  the  elements,  the 
structural  stiffness  matrix  becomes,  for  N  elements 

N  m 

[K]  =  Z  2 iff  [B]  [D][B]rdrdz  (A-3! 

n=l 

Since  [B]  contains  terms  dependent  on  the  coordinates, 
straightforward  integration  may  not  be  possible.  An  approxi¬ 
mate  solution  can  always  be  obtained  by  evaluating  [B]  at  a 
centroid  defined  by 

r  =  (r^  +  r2  +  r3)/3  (A-4a! 


z 


(z,  +  +  zJ/3 


(A- 4b 


APPENDIX  B 


This  section  contains  experimental  data  for  the  Centrum 
tests.  It  includes  Calibration  Data,  Experimental  Curves, 
Extracted  Data,  and  Internal  and  External  Measurements. 


TABLE  B-l 


Calibration  Data* 


Specimen 

I.D. 

Load 

(lb) 

Calibration  Factor 
(10“4  in/in) 

Zero 

(in) 

KAZ  7  (LI) 

15 

0.9366 

-8 

KAZ  7  (L2) 

15 

0.9198 

+4 

KAZ  10 (LI) 

15 

0.9901 

+7 

KAZ  10 (L2) 

15 

0.9671 

-3 

KAZ  30 (LI) 

15 

0.9398 

+10 

KAZ  30 (L2) 

15 

1.0024 

+8 

KAZ  3 (LI) 

30 

1.0378 

+7 

KAZ  3  (L2) 

30 

0.9634 

+9 

KAZ  6  (LI) 

30 

0.9679 

-4 

KAZ  6(L2) 

30 

1.0229 

0 

RHES  17 (LI) 

30 

0.9992 

+14 

RHES  17 (L2) 

30 

0.9785 

+32 

*  Displacement  (in)  =  Calibration  Factor  x  (Digitized  Value  - 
Zero) 
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Midplane  displacement  is  half  of  specimen  displacement. 

Data  from  creep  plots.  Specimens  are  identified  by  primate  name  and  location  (in  parentheses) . 
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Fig.  B-5.  Test  curve,  specimen  KAZ  6 (L2) 
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Test  curve,  specimen  RHESUS  17 (Ll) 


APPENDIX  C 


Example  of  Calibration  Procedure 

This  appendix  presents  an  example  of  the  procedure  used 
to  determine  the  calibration  factor  described  in  section 
3.1.2  of  the  text. 

Consider  the  following  successive  displacement  readings 
taken  at  the  data  terminal  for  the  indicated  cylinders: 

Cylinder  Used  Reading 

Long  -2712 

Short  -1413 

Long  -2699 

The  average  reading  using  the  long  cylinder  is  -2705.5.  The 
measured  difference  between  this  average  reading  and  the 
reading  from  the  short  cylinder  is  therefore  -1292.5.  Finally 
the  calibration  factor  is 

..  .  0.1250  ..  .. 

Calibration  Factor  =  :  7  (3-1) 

[Measured  difference!  10 

_  0.1250 

1-1292. 5|  10"4 

Calibration  Factor  =  0.9671 


APPENDIX  D 


Method  for  Predicting  q-^ 

This  appendix  describes  a  method  for  graphically  pre¬ 
dicting  an  approximate  value  for  the  material  parameter 
described  in  this  thesis.  It  is  based  on  a  knowledge  of  the 
axial  displacement  -  time  behavior  immediately  after  creep 
strain  begins. 

Consider  the  axial  displacement  -  time  curve  illustrated 
in  Fig.  D-l,  where  uZq  represents  the  initial  elastic  dis¬ 
placement  at  t  =  tQ/  and  uz^  represents  the  displacement  at 
t  =  t^,  a  short  time.  At,  later.  An  approximate  expression 
for  the  initial  displacement  rate  uQ  is  the  slope,  a,  of  a 
straight  line  between  the  two  points,  where 


a  = 


(D-l) 


An  approximate  expression  for  creep  strain  rate,  ec,  can  thus 
be  obtained  as 


(D-2) 


where  Lq  is  an  average  axial  length  immediately  after  elastic 
deformation  has  occurred  and  before  any  creep  strain  has 
accumulated.  Using  the  fact  that  is  related  to  the  slope 
of  the  stress  versus  strain  rate  curve  for  cc  =  0  (see  the 
main  text  for  a  discussion) ,  and  using  an  average  axial  stress 


°AVG'  over  the  loadad  surface  at  t  =  t  ,  we  can  approximate 

qx  by 

q,  =  °AVG  (D-3) 

1  «C 

Estimation  will  improve  as  At  is  made  smaller,  since  the 
true  value  of  will  then  be  approached  in  the  limit  as  At 
approaches  zero.  This  approximate  value  can  then  be  used  as 
a  good  first  guess  in  which  to  begin  the  analytical  technique 
described  in  this  thesis.  The  author  found  it  useful,  per¬ 
mitting  fewer  computer  runs  to  converge  to  a  final  value. 
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